%-------------------------------------------------------------------------
% Sub routine: Estimation and identification
%
% Estimation Codes for State-Level Results in  
% Ben-David, Itzhak, Sebastian Weber, and Pascal Towbin 
% "Inferring Expectations from Observables: Evidence from the Housing Market" 
% Review of Economics and Statistics
% ------------------------------------------------------------------------ 

%% Start loop over states
for is=1:51
if is==37 % Drop state with msising data
else
    
% state identifiers
ct=is;
state=State(t*ct,1);
state_vec(ct,1)=state;

%% Find start and end observation and construct estimation data
construct_y_sta

%% Estitmation 
tau=0;   % uninformative prior
[gamma,Su,resid,Par] = BVARgibbs_parameters(y,p,tau); 

%% Identification
% initiate
tel=0;
tel2=0;
rot=0;
nonstab_tot=0;

while tel<n_rot_acc  

jj=1;
    
[choleskyimp,A(:,:,ct) ] = IRF_BVARGibbs_2(gamma,Su,p,n_var,1,n_step); %% Sampling 
gamma_g(:,:,ct)=gamma;
  
rot=rot+1;
    
           
%Draws Rotation matrix
W = normrnd(0,1,n_var,n_var);
Q = getqr(W);

      for tt=1:n_step
            irfcand(:,:,tt)=choleskyimp(:,:,tt,jj)*Q;  
      end
      
irfcand(n_pr_cum,:,:)=cumsum(irfcand(n_price_n,:,:),3);

J = zeros(n_var,n_var*p); J(:,1:n_var) = eye(n_var);

% impose sign restrictions
        for shock=1:n_var
            
            % (1) Supply shock
         
               test_su= [checksign3(irfcand(n_pr_cum ,shock,:),LL,KK,'+')  
                checksign3(irfcand(n_total_vac, shock,:),LL,KK,'-')   
                checksign3(irfcand(n_inv_l ,shock,:),LL,KK,'-')
                ];
          
             if min(test_su)==1
                 shockcheck(shock,1)=1;
             elseif max(test_su)==-1;
             shockcheck(shock,1)=1;
             irfcand(:,shock,:)=-irfcand(:,shock,:);
             else
                shockcheck(shock,1)=0;
             end
             
             % (2) Housing consumption shock
         
              test_de=[  checksign3(irfcand(n_pr_cum, shock,:),LL,KK,'+')
                checksign3(irfcand(n_inv_l ,shock,:),LL,KK,'+')
                checksign3(irfcand(n_total_vac, shock,:),LL,KK,'-')  
                checksign3(irfcand(n_mrate, shock,:),LL,KK,'+') 
                      ];
      
             if min(test_de)==1
                 shockcheck(shock,2)=1;
             elseif max(test_de)==-1;
             shockcheck(shock,2)=1;
             irfcand(:,shock,:)=-irfcand(:,shock,:);
             else
                shockcheck(shock,2)=0;
            end
             
            % (3) Mortgage rate shock

            test_int=[
                checksign3(irfcand(n_pr_cum,shock,:),LL,KK,'+')
                checksign3(irfcand(n_inv_l ,shock,:),LL,KK,'+') 
                 checksign3(irfcand(n_mrate ,shock,:),LL,KK,'-')   
               ];
           
          if min(test_int)==1
                 shockcheck(shock,3)=1;
             elseif max(test_int)==-1;
             shockcheck(shock,3)=1;
             irfcand(:,shock,:)=-irfcand(:,shock,:);
             else
                shockcheck(shock,3)=0;
          end
            
          
            % (4) Expectation shock

            test_sp=[
                checksign3(irfcand(n_pr_cum,shock,:),LL,KK,'+')          
                checksign3(irfcand(n_inv_l ,shock,:),LL,KK,'+')     
                checksign3(irfcand(n_total_vac, shock,:),LL,KK,'+')  
                checksign3(irfcand(n_mrate ,shock,:),LL,KK,'+')      
                ];
            
            if min(test_sp)==1
                 shockcheck(shock,4)=1;
             elseif max(test_sp)==-1;
             shockcheck(shock,4)=1;
             irfcand(:,shock,:)=-irfcand(:,shock,:);
             else
                shockcheck(shock,4)=0;
            end
            

        end

        controlvec= sum(shockcheck,1);
              
        if controlvec==ones(1,size(controlvec,2));
             tel=tel+1;
             if tel/10==fix(tel/10)
             disp([num2str(tel) ' accepted draws out of ' num2str(rot) ' draws']);
             end
            [irfstr(:,:,:,tel,ct), vardecomp(:,:,:,tel,ct)] = irfselec(shockcheck,irfcand);
            BB(:,:,tel)=irfcand(1:n_var,1:n_var,1);
           errstr=BB(:,:,tel)\resid';
            strerr(:,:,tel,ct) = ord_strerr(shockcheck, errstr);
            
           
        end
    
   clear irfcand
   rot=rot+1;
    end

end
irfstr_nn=irfstr;
vardecomp_mean=mean(vardecomp(:,:,:,:),4);
end
no_str=size(shockcheck,2);

% %% Transform into level  and save IRFs and errors
irfstr_nn(cum,:,:,:,:)=cumsum(irfstr(cum,:,:,:,:),3);
irfstr_g4=irfstr;
strerr_save=strerr;
